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Preface 


A research effort for surface radiation budget (SRB) studies was initiated at the NASA Langley 
Research Center (LaRC) in the mid-1980’s with the goal of developing parameterized, fast radiative 
transfer algorithms for deriving shortwave (SW) and longwave (LW) components of SRB on a global 
scale using meteorological products from operational satellites. During the ensuing years, a group of 
scientists at the LaRC developed SW and LW algorithms and derived global fields of all components of 
SRB. The product of this effort was the first long-term global dataset of SRB components based on 
satellite data. This dataset has since been published and made available to the climate science 
community. 

The SW algorithm for the above effort was developed by W. F. Staylor (with help from W. L. 
Darnell and the first author of this report), and came to be known as the Staylor algorithm. Staylor 
accounted for extinction of solar radiation in the atmosphere by adopting existing parameterizations from 
the literature for some processes, modifying existing ones for others, and developing new ones for still 
others. The early versions of this algorithm used meteorological products from the TIROS Operational 
Vertical Sounder (TOVS) and the Advanced Very High Resolution Radiaometer (AVHRR) flown aboard 
NOAA’s operational satellites, and were fine-tuned using surface insolation measurements from a number 
of sites in the U. S. Later versions were adapted to use meteorological data and other inputs from the 
International Satellite Cloud Climatology Project (ISCCP) and the Earth Radiation Budget Experiment 
(ERBE), and were extensively validated with surface measurements from the Global Energy Balance 
Archive (GEBA) and other sources. 

Despite its simplicity and computational efficiency, this SW algorithm performed remarkably 
well when compared with much more detailed and computationally slower algorithms. It compared 
favorably in the InterComparison of Radiation Codes in Climate Models (ICRCCM) and other validation 
exercises sponsored by the World Climate Research Program (WCRP). This led to its selection by the 
WCRP/SRB project in the early 1990’s as one of the two algorithms for producing global insolation 
datasets. It was also selected by the Global Energy and Water-cycle Experiment (GEWEX)/SRB 
Workshop in 1993 as a quality-check algorithm for monitoring the performance of the primary algorithm 
chosen by the GEWEX/SRB project. Currently, it is being tested by the NASA/Clouds and the Earth’s 
Radiant Energy System (CERES) project along with other surface SW algorithms. There is considerable 
interest on the part of the CERES project hi using this algorithm for deriving surface SW products. 

Along with the strengths mentioned above, this algorithm suffered from a number of serious 
weaknesses. Generally, the sources of information used were not well documented. Parameterizations 
supposedly taken from the literature were difficult to trace to their referenced sources. The assumptions 
and approximations made in the algorithm were not adequately explained or justified. Seemingly 
complicated formulas developed for the algorithm were not supported by a satisfactory analytical 
framework. These shortcomings often resulted in less than enthusiastic acceptance and use of this 
algorithm by the scientific community. 

Despite the noted shortcomings, superior performance of this algorithm sustained the interest of 
the GEWEX/SRB and CERES projects in its continued use. To meet the needs of these projects, another 
effort was undertaken at LaRC to fully document the scientific basis of this algorithm. Efforts were made 
to relate the information taken from the literature to its original sources, clarify and justify the 
assumptions and approximations used, and establish the analytical framework for the formulas developed 
and used. Changes were made when some components of the algorithm could not be justified. Also, a 
fully documented FORTRAN-90/95 code was developed for implementing this algorithm using newer 
sources of meteorological inputs and providing results at higher spatial resolution. The present report 
represents the first step of the latter effort. Having undergone large-scale changes during this process, and 
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with many more changes on the way to meet the future needs of the new projects, this algorithm needed a 
new identity. It has been renamed the Langley Parameterized Shortwave Algorithm (LPSA). 
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1. Introduction and background 

The surface radiation budget (SRB) is a major component of the energy exchange between the 
atmosphere and the land/ocean surface and thus exercises a profound influence on many weather and 
climate processes (Ramanathan 1986). Developing a long-term climatology of SRB on a global scale was 
recognized as an essential prerequisite for a number of World Climate Research Program (WCRP) 
research projects (Suttles and Oliring 1986). In response, the WCRP established the SRB Climatology 
Project (WCRP-69 1992) to facilitate and monitor the development of a reliable long-term climatology of 
SRB. Recognition of the scientific potential of SRB also led to the creation of a SRB program at the 
NASA Langley Research Center (LaRC) starting in mid-1980’s. The main objective of the LaRC 
program was to develop efficient algorithms for producing SRB parameters on a global scale, preferably 
using input data from operational satellite sources. 

The algorithm for deriving surface SW fluxes for the LaRC program was developed by W. F. 
Staylor and colleagues (Staylor et al. 1983; Darnell et al. 1988, 1992), and became known as the Staylor 
algorithm. Meteorological products from NOAA’s polar-orbiting satellites were chosen as inputs for the 
algorithm because of their global coverage. The first version of this algorithm used water vapor and 
ozone abundances from TO VS data products. Planetary albedos derived from AVHRR data were used as 
proxy for cloud parameters. The second version was tailored to take advantage of the Cl datasets as 
those became available from ISCCP (Rossow and Schiffer 1991). Cloud products derived from the 
ISCCP network of geostationary satellites, intercalibrated with AVHRR data taken from one or two polar 
platforms, resulted in a much improved SRB dataset from the second version. A global SRB climatology 
derived using the second version and the entire 8-year record of ISCCP-C1 data (July 1983 to June 1991) 
has been published recently (Gupta et al. 1999). Also, this dataset has been made available online to the 
science community worldwide from the LaRC Atmospheric Sciences Data Center 
(http://eosweb.larc.nasa.gov). 

The purpose of this report is to present a new version of the Staylor algorithm with detailed 
description and documentation of its methodology. Also, a FORTRAN-90/95 code has been developed 
for the implementation of the new version. Section 2 identifies the limitations of the older versions and 
makes a case for restructuring the algorithm for newer applications. Section 3 examines the simple basis 
of the algorithm and links that basis, wherever possible, to its original sources in the literature. Section 4 
describes the sources of input data used in the past, in use at present, and being considered for the future. 
A complete list of surface SW parameters and sample results derived from the new version are presented 
in Sec. 5. A short list of planned enhancements is presented in Sec. 6. 

2. Case for restructuring 

The Staylor algorithm is simple, computationally efficient, and above all, provides good results. 
Its potential remains largely unrealized, however, because it has not been accepted widely by the science 
community. It is not difficult to identify the reasons for this non-acceptance, the foremost being the lack 
of detailed documentation of the methodology. In some cases, material referenced from literature sources 
was not traceable to those sources, or the assumptions and approximations used were not sufficiently 
explained or justified. In other cases, complex formulas developed for the algorithm were not supported 
by an adequate analytical framework. It became obvious that such conditions had to be remedied before 
this algorithm and its results would gain wider acceptance. 

There are a number of other reasons which motivated the development of the new version of this 
algorithm. Based on its superior performance, the WCRP/SRB (same as the GEWEX/SRB) project had 
designated it as the quality-check algorithm for monitoring the performance of the primary SW algorithm 
developed by Pinker and Laszlo (1992a). The GEWEX/SRB project is required to produce surface fluxes 
on a 1° x 1° spatial resolution using meteorological inputs from data assimilation models. Cloud 
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parameters for the GEWEX/SRB project are being produced on a 1° x 1° spatial grid using pixel-level 
(DX) data from ISCCP (Stackhouse et al. 2000). The CERES project at NASA/LaRC currently produces 
net SW fluxes at the surface using satellite-derived top-of-atmosphere (TOA) broadband SW fluxes with 
the TOA-to-surface transfer algorithm developed by Li et al. (1993). Validation of net SW fluxes is 
difficult to accomplish because the most commonly measured surface SW parameter is insolation. 
Algorithms for converting net SW flux to insolation are of limited use because they require prior 
knowledge of surface albedo fields. As a result, there is interest in the CERES project in using this 
algorithm for producing surface insolation independently of the Li et al. algorithm. Continued interest by 
the GEWEX/SRB and CERES projects in the use of this algorithm necessitated the development of a 
well-documented and scientifically justifiable version. This new version has been named the Langley 
Parameterized Shortwave Algorithm (LPSA). 

3. The algorithm 


On a very basic level, surface insolation (downward SW flux at the surface), F SD , in the LPSA is 
computed as 


F = F T T 

r SD r TOA 1 A 1 C 


( 1 ) 


where F TOA is the corresponding insolation at the TOA, T A is the transmittance of the clear atmosphere, 
and T c is the transmittance of clouds. All quantities in Eq. (1) refer to the broadband SW region, 
approximately from 0.3 to 5.0 pm. The temporal resolution of the current version of LPSA is daily, i.e. 
all parameters are computed on a daily average basis. 


3.1. TOA insolation 


F toa is computed as (Peixoto and Oort 1993) 

f toa = S(d m /d) 2 cos Z, 


( 2 ) 


where: 

S - solar constant, 

d - instantaneous Sun-Earth distance, 

d m - mean Sun-Earth distance, also called the astronomical unit, and 
Z - solar zenith angle. 

The value of cos Z for any location and time is computed as 

cosZ = sin0 sind + cos0 cosd cos/?, (3) 

where <p is the latitude of the location, <5 is the solar declination at the time, and h is the hour angle (in 
radians) from the local meridian. For a more detailed definition of /?, the reader is referred to Sellers 
(1965). A few simple relationships that follow from Eq. (3), and will be of interest later, are presented 
below. At sunrise and sunset, Z = jt /2, and h = H (defined as half day length), so that 

sin 0 sin <5 + cos <p cos <5 cos// = 0, (4) 


which gives 


4 



cos H 


- tan U tan/7, 


(5) 


or 

H = cos -1 (- tan{7tan/7)- (6) 

The eccentricity correction parameter, (djdf, and the solar declination, [J are both functions of 
time. Strictly speaking, Eqs. (2) and (3) used for computing F TOA are applicable only on an instantaneous 
basis. In the current work, however, simple expressions presented by Iqbal (1983) for computing daily 
average values of (djdf and 0 w'ere used to match the daily temporal resolution of the algorithm. A brief 
discussion of those expressions and the maximum errors of daily average values computed with them is 
presented in Appendix A. The maximum error incurred in computing F TOA using daily averages instead of 
instantaneous values of (djdf and 0 was estimated to be less than 0.01%. 

Daily total insolation at the TOA, FD TOA , can be computed by integrating Eq. (2) from sunrise to 
sunset as 

sunset 

FD TOA = S { d m/ d f f COS Z dt . (7) 

sunrise 


When time is measured in hours, 


, 12 „ 
dt = — dh, 
0 


( 8 ) 


and 


FDtoa ~ u S{ d m/ d ) or sin/7sin/7d7; + cos 0 cos 0 cos h dlij , (9) 


Integrating Eq. (9) gives 


0/1 

FDtqa = — S (d m ! d)~ (H sin/7sin/7+ cos/Jcos/Jsin//), 


( 10 ) 


and hourly average TOA insolation, FDA toa as 


FDA toa = — S(d m /d)" (//sin/7 sin/7 + cos/7 cos/7 sin//), 


0 


( 11 ) 


averaged over a 24-hour day. 


Staylor computed daily average TOA insolation in tenns of a parameter D, defined as 


D 


D 


F cos" 1 (-FIG) + G Vl -(FIG) 1 


( 12 ) 


where F = sin /7 sin Q, and G = cos jj cos [j Equation (12) is the same as Eq. (11) without the Sun -Earth 
distance corrected solar constant, S (djdf. Staylor named the parameter D as the daily mean vertical Sun 
fraction. D is the ratio of the actual daily TOA insolation to that if the Sun was overhead for the entire 24 
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hours. D can also be viewed as the value of cos Z averaged over the 24-hour day. Staylor also defined 
another parameter, u, as the average of cos Z from sunrise to sunset (Staylor and Wilber 1990), in the 
fonn 


u = F + G^(G[]F)/2G. (13) 

The authors were not able to derive an expression for u in the form shown in Eq. (13). It was decided, 
therefore, to derive an expression for u from first principles as 


j sunset 

u = ricos Z dt , 

day length I— I 

sunrise 


(14) 


where ‘ day length ’ represents the time from sunrise to sunset in hours. Equation (14) can also be written 
in terms of F, G, and h as 


u = -jj; Q ( F + G cos h)dh = F + t 


or 

Fcos ni (UF/G) + G^ IU(F/G) 2 
U cos D1 (UF/G) 

The values of u computed with Eqs. (13) and (16) were found to be significantly different. The highest 
value of u which occurs at the equator during the equinoxes was found to be 0.707 for Eq. (13) and 0.637 
for Eq. (16). The parameter u was extensively used by Staylor for deriving approximate values of a 
number of input variables (e.g., see Eqs. 39, 42, 43, 46, and 49) which are known to be dependent on Z 
(or cos Z). These input variables were essential for computing insolation but were not always available 
for all conditions from their regular sources. This strategy proved valuable in eliminating large gaps in 
the computed flux fields when those input variables were not available. 

Equations (12) and (16) for D and u respectively, apply over most of the globe where there is 
sunrise and sunset during the course of a day, and —1 < F/G < 1. Over polar regions of the summer 
hemisphere, there is no sunset, F/G > 1, and 

cos d1 (DF/G) = H = 0 - (17) 

Tire tenn G ^ 1 0 (F/G) 2 in Eqs. (12) and (16) becomes undefined. At those latitudes, 

D = u = F. (18) 


(15) 


(16) 


6 



Over polar regions of the winter hemisphere, there is no sunrise, FIG < -1 , both D and u are undefined 
and are set equal to zero. 


3.2. Clear-sky transmittance 

Clear-sky transmittance, T A , was computed as 

T a = (1 + B)e ar z, (19) 

where B is a backscatter tenn (defined later in Eq. 33), t z is the broadband extinction optical depth at 
solar zenith angle Z and accounts for all absorption and scattering processes in the clear atmosphere. 
Staylor analyzed clear-sky insolation measurements from several sites within the United States to 
empirically derive the dependence of r z on Z in the form 

r z = t 0 (sec Z) N , (20) 

where r 0 is the broadband optical depth for overhead Sun (Darnell et al. 1992). In first two versions of 
the algorithm, an expression for the exponent N was derived empirically in the form 

N = 1.10 2.0 r 0 . (21) 

Another method, to be presented later in this section, is being used for deriving N in the current version of 
the algorithm. 

The vertical optical depth (for overhead Sun), r 0 , was computed as 

r 0 = □ In (1 0« o ), (22) 

where a 0 is an effective vertical attenuation factor for broadband radiation. Functionally, a 0 is the 
equivalent of an absorptance, and may be called an ‘extinctance.’ Staylor computed a 0 in terms of 
effective attenuation factors for the many absorbing/scattering constituents of the clear atmosphere as 


«n 


= a 


H-,0 


+ a n + a, 


co , 


+ a n + a 


Ray 


+ a 


Aer ’ 


(23) 


where the first four terms on the right represent attenuation by the respective absorbing gases, a Ray 
represents attenuation due to Rayleigh scattering, and a Aer , the attenuation due to aerosols. Staylor’ s 
choice of this formulation was guided by the fact that all attenuation factors except a Aer were already 
available in the literature. Note that the addition of individual attenuation terms as shown in Eq. (23) is 
justifiable only if the processes represented by them occur independently of one another. In simple terms, 
it means that either these processes occur in different regions of the spectrum, or in different altitude 
ranges of the atmosphere. Among the processes dominating the extinction of solar radiation in the 
Earth’s atmosphere, absorption by ozone takes place in the ultraviolet and visible regions and primarily in 
the stratosphere. Rayleigh scattering also occurs in the ultraviolet and the visible but mostly in the 
troposphere. Water vapor abosorption takes place in the near infrared and in the troposphere. These 
processes may be considered independent of one another in accordance with the above requirement. Tire 
remaining extinction processes are much less significant. 
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Staylor used the following formulas to represent attenuation factors for the various processes: 


0h 2 o 

= 0.100 (U H2 q) 0 - 21 , 

(24) 

0o 3 

= 0.037 (U O} ) 0A3 , 

(25) 

1 

co 2 

0.006 (P s 350 ) 029 , 
5 300 

(26) 

0o 2 

= 0.0075 ( P s ) 0 87 , 

(27) 

l-l Ray 

= 0.035 (P s ) °' 67 , 

(28) 


where U H20 represents the column water vapor amount (in precipitable cm), U 03 , the column ozone 
amount (in cm-atm), and P s , the surface pressure (in atm). These formulas have changed only slightly 
between evolving versions of the algorithm. The attenuation factors in Eqs. (26) — (28) are tied to the 
surface pressure which acts as a surrogate for column abundances of the uniformly mixed gases. The 
fraction (350/300) in Eq. (26) represents a scaling of C0 2 attenuation to the present-day mixing ratio of 
350 ppm from the 300 ppm most likely used in the original formula. Equation (27) above was taken from 
the first version of the algorithm (Darnell et al. 1988) even though 0.002 (in place of 0.0075) appeared in 
later versions. The form of Eq. (27) given above matches exactly with the formula for 0 2 absorption 
given by Hoyt (1978). Staylor attributed attenuation factors of Eqs. (24), (25), and (28) to Lacis and 
Hansen (1974), and those of Eqs. (26) and (27) to Yamamoto (1962). An examination of the above 
papers shows that formulas in the fonns as in Eqs. (24) - (28) are not presented in those papers. A closer 
examination of Eqs. (24) - (28), and their relationships and comparisons with the material contained in 
those papers is presented in Appendix B. 

The fonnula used for U Aer in the current version of the algorithm has the form 

0 Aer = O Aer ^UoJ 0 ) + ^O Aer (o 0 (lUg), (29) 

where /X,, co o> and g are the broadband aerosol optical depth, single scattering albedo, and asymmetry 
parameter, respectively. Tire first term on the right represents attenuation by aerosol absorption, and the 
second tenn, by aerosol backscattering (Wiscombe and Grams 1976). Forward scattering by aerosols 
does not cause attenuation of the downward radiation stream and, therefore, is not included in Eq. (29). It 
is noted that Staylor did not use u> 0 in the second term assuming a> 0 to be always close to unity. The 
lowest value of a) 0 used by Staylor (for continental aerosol, see Appendix C) was 0.90. Information on 
the geographical distribution and radiative properties of aerosols used in the current work is presented in 
Appendix C. Also, for the sake of completeness, it needs to be mentioned that in earlier versions (Darnell 
et al. 1992), the fonnula for jJ Aer had the form 

0 Aer ~ 007 + 0.009 (U Hi0 ), (30) 

which tied [ ] Aer to the moisture content of the atmosphere. 

The value of the exponent N of Eq. (20) for the current version was derived using computed at 
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two values of Z. Z = 0° and 70.5° were chosen to cover a wide range. Corresponding values of [J are 

denoted as jj, and [ ~} 0 . [] 0 was computed from Eq. (23), and Q> therafter from Eq. (22). Values of jj 70 and 

[Jo were computed by substituting 3 times the constituent abundances in Eqs. (24) — (29) and then 
applying Eqs. (23) and (22). Note that for Z = 70.5°, sec Z, which represents the relative air mass is equal 
to 3. Finally, from Eq. (20) 

0o = a (sec Z) N = Oo( 3) n , (31) 

which gives 

N = — { ln(/X,) D l n (Do) } > (32) 

m3 

providing an average value of N for the entire range of Z from 0° to 70.5°. 

The term B introduced in Eq. (19), represents radiation backscattered by the atmosphere after 
being reflected upward from the surface. It represents an enhancement of the downward radiation stream 
and currently has the form 

B = 0.065P s A s + 2A s n Aer D 0 ^Ug), (33) 

where A s is the surface albedo. The first tenn on the right represents the Rayleigh backscattering of the 
surface reflected (diffuse) radiation. The value of the coefficient in the first term (0.065) is in good 
agreement with the estimate of spherical albedo of a Rayleigh atmosphere illuminated from below (Lacis 
and Hansen 1974). The second term represents the surface reflected radiation backscattered by aerosols. 
In this term also, Staylor did not use f] 0 for the same reason as stated in connection with the second tenn 
of Eq. (29). Note that the magnitude of this enhancement in Eq. (33) is four times larger than the 
attenuation by a similar process represented by the second term in Eq. (29). This is a result of two causes: 
i) the surface reflected radiation is completely diffuse and the scattering optical depth for diffuse radiation 
is approximately double that for the direct radiation (Wiscombe and Grams 1976), and ii) the reflected 
radiation makes two passes through the atmosphere doubling the scattering optical depth again. Note 
further that the aerosol term was not included in the expression for B in the first two versions of the 
algorithm. In the second version (Darnell et al. 1992), only the first tenn on the right in Eq. (33) was used 
to represent B. Prior to that, in the first version (Darnell et al. 1988), the backscattering enhancement was 
included as a negative term (-0.065 P s A s ) in the expression for jj 0 in Eq. (23). 

3.3. Cloud transmittance 

Cloud transmittance, T c , for a grid box was computed in tenns of visible reflectances as 

T r = 0.05 + 0.95 <7?ovcD7? " ,p " ) , (34) 

(*«□**) 

where R oyc , R ctn and R meas represent daily average values of overcast, clear, and measured reflectances 
respectively corresponding to an overhead Sun. Values of R meas were computed by averaging measured 
daytime instantaneous (3-hourly) reflectances, weighted by pi 0 , the instantaneous value of cos Z. Equation 
(34) is based on standard threshold methods used for cloud parameter determination (e.g., Moser and 
Raschke 1984) and a recognition of the observational fact that even for the thickest clouds, T c is not 
reduced to zero. 
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Since completely overcast conditions do not occur every day, R orc was computed using a 
regression relation of the form 


R 


Q 


+ 


C, 


M/A M if, 


M /'„ 

( M +/0 


(35) 


where p is the cosine of the instantaneous view zenith angle. Coefficients C. t and C 2 were determined off- 
line, separately for each ISCCP satellite for every month, by linear regression between R ovc fj,fj, B and 
\uuj(u+u 0 )\ > . The theoretical basis for Eq. (35) can be found in Staylor (1985). Reflectances used in the 
above regression were sampled from the monthly ensemble of overcast reflectances on the basis of cloud 
optical depths. Only reflectances corresponding to the highest 10 — 20% range of cloud optical depth 
values were selected. 


At least two different methods were used for computing the values of R clr depending on the 
underlying surface type. For ice-free ocean surface, R clr was computed from 

R lr = C 3 + C 4 (ppj 00 ' 75 . (36) 

Information on the theoretical basis of Eq. (36) was not available. Coefficients C s and C 4 were also 
determined for each ISCCP satellite for every month by linear regression between sampled daily averages 
of R ch . and (pp„) 0 75 for clear-sky grid boxes over ice-free ocean. An alternative method was used for all 
surface types other than ice-free ocean. With this method, R ch . for each grid box was computed by 
averaging all available daytime values of clear-sky visible reflectance. The authors note that even though 
different methods of computing R ch . were suggested in Darnell et al. (1992) for different surface types 
(e.g., land, snow-covered land), only the method described in this paragraph was implemented in the 
code. Details of the sampling procedures used for both regression analyses (Eqs. 35 and 36) are presented 
in Appendix D. 

It should be noted that the use of Eq. (34) was not found to be appropriate for all grid boxes and 
all days. Under certain conditions, other methods for computing T c had to be used. A discussion of the 
conditions under which alternative computation of T c was necessary and the equations used for those 
computations is also presented in Appendix D. 

3.4. Surface albedo 


Surface albedo, A s , is an important SRB parameter as the primary determinant of net or absorbed 
SW radiation, FsNi which is computed as 

F sn = F SD (1[]A S ). (37) 

On a secondary level, surface albedo also affects the downward SW radiation through the backscattered 
radiation term represented by Eq. (33). The all-sky surface albedo used in Eq. (37) was computed by 
Staylor as 


A? - Asovc + ( Asc/r D Asovc ) ’ (38) 

where A Sclr and A Sovc represent surface albedos for clear-sky and overcast conditions respectively. A Sclr and 
A Smc may be substantially different because of the differences between illumination geometry under clear- 
sky and overcast conditions. Staylor obtained A Sclr and A Soyc from different sources for different surface 
types as described below. 
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Ice-free oceans: For ice-free oceans, Staylor (Darnell et al. 1992) computed surface albedos as 


A Sclr = 0.039 !u, 


( 39 ) 


and 


^ Save = 0-065. (40) 

Even though Darnell et al. (1992) refer to Payne (1972) and Kondratyev (1973) as the sources of Eq. (39), 
an equation of this form is not found in those documents. Closer examination of those documents shows, 
however, that the results reported therein agree with those represented by Eqs. (39) and (40). The authors 
believe that Staylor developed Eq. (39) in its simple form using results reported in Payne (1972) and 
Kondratyev (1973). Darnell et al. (1992) also showed that Eq. (40) follows from Eq. (39) for u = 0.60, 
which in turn, corresponds to Z = 53°. Note that 53° represents a good estimate of an effective value of Z 
for diffuse radiation. 

Other snow/ice-free surfaces: For all snow/ice-free surfaces other than ice-free ocean, Staylor 
used surface albedos derived from monthly-average clear-sky TOA albedos obtained from the Earth 
Radiation Budget Experiment (ERBE; Barkstrom et al. 1989). ERBE-based surface albedos were used 
whenever and wherever they were available, and were derived as described below. Linear relationships 
between clear-sky TOA albedo, A Tck , and corresponding surface albedo have been developed over the 
years in the simple form 

Kir = a + b A Sclr , (41) 

(e.g., Chen and Ohring 1984; Koepke and Kriebel 1987), generally on an instantaneous basis. The 
constants a and b represent the effect of the intervening atmosphere and are functions of the atmospheric 
properties. Staylor (see Staylor and Wilber 1990) adapted Koepke and Kriebel’s version of Eq. (41) for 
daily average values of A Tclr and A Sdr . This version includes the effects of Rayleigh scattering, water 
vapor, ozone, and aerosols built into the constants a and b. For the daily average form of Eq. (41), 
Staylor and Wilber (1990) represented the above constants as 

a = 0.25 P s / (1 + 5 m ), (42) 


and 


b = 1 □ a □ 0.04 | l6U 0i /(1 + 5 m )| 06 □ 0.12(U H2 q ) 025 
□ 2.40 Aer ^UD 0 )u° A U n Aer 1(2 + 15 u L5 ), 


and computed the value of A Sclr as 


A sdr = (ArcHa)/ b. (44) 

Note that expressions for a and b as given in Eqs. (42) and (43) were not found in either Chen and Ohring 
(1984) or Koepke and Kriebel (1987). Also, in earlier versions of the algorithm (see Darnell et al. 1992), 
a simpler relationship of the fonn 


Asdr ~ 1 -5 A lch . G 0.07, 


(45) 
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was used which was based primarily on the work of Chen and Ohring (1984). 

Monthly average clear-sky TO A albedos ( A Tch . ) from ERBE were available for a period of 51 
months from March 1985 to May 1989. Corresponding monthly average values of A Sclr for these months 
were computed from Eq. (44) and used as such for every day of the respective months. The possibility of 
using ERBE-derived values of A Sclr for months outside the ERBE period was also examined. To that end, 
interannual variability of ERBE-derived A Sch . for each month over the available years was analyzed. This 
analysis showed interannual variability of up to 10% over high and mid latitudes of the Northern 
Hemisphere (NH), and in the 1 — 2% range over lower latitudes and most of the Southern Hemisphere 
(SH). Staylor, therefore, decided to use multi-year averages of A Sch . for each month from the ERBE 
period, for corresponding months outside the ERBE period. This practice is being continued until better 
surface albedo datasets become available. For all snow/ice-free regions where ERBE-derived A Scll . was 
used, A Smc was derived from 

As™ = 1-1 A Sclr u 0 - 2 . (46) 

Surfaces affected by snowlice: Staylor also made use of ERBE-derived surface albedos for 
regions which were affected by snow/ice. For such regions, however, ERBE-derived values were 
modified to account for the presence of snow/ice. When ERBE-derived values were not available for 
some snow/ice affected regions, Staylor devised other relationships for deriving surface albedos using 
snow/ice fractional cover for that region and completed the flux calculation. The forms of those 
relationships and the conditions under which they were used are discussed in Appendix E. 

It is important to note here that Staylor’ s choice of Eqs. (39) and (40) over oceans, and ERBE- 
based surface albedos over other regions helped overcome a serious difficulty in deriving broadband SW 
fluxes. These albedos were already broadband. Attempts to derive broadband surface albedos from 
ISCCP visible radiances were plagued with large uncertainties involved in the narrowband-to-broadband 
conversion. 

3.5. Direct, diffuse, and PAR 

Direct and diffuse broadband fluxes, and PAR (photosynthetically active radiation, between 0.4 
and 0.7 pm) are components of global insolation which are important for a variety of applications. 
Staylor devised simple formulas for partitioning global flux into direct and diffuse components based 
primarily on cloud transmittance, T c . Specifically, the partitioning formulas depended on whether or not 
a value for T c for the grid box was available. When T c was available, and it was > 0.35, the direct flux at 
the surface, F sdin was computed as 

F sdir = F S d(TcU 0-35), (47) 

and the diffuse flux at the surface, F sdi f, as 

Fsdif = Fsd ( 1 -35 UT C ). (48) 

Together, these equations mean that for clear skies (T c = 1.0), direct flux is 65%> of the global flux, and 
the remaining 35% is diffuse flux. Further, when T c [] 0.35 (generally dense cloudiness), the direct flux is 
reduced to zero and the entire global flux is in the diffuse form. Staylor almost always devised an 
approximate method when a variable required to complete the calculation was not available from its 
standard sources. In keeping with that approach, Staylor adopted the 65/35 partitioning (the same as for 
clear skies) when a value of T c could not be determined from the usual methods. Exact details on how 
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Staylor derived Eqs. (47) and (48) were not available, but the authors believe that these equations were 
developed by fitting curves to the results obtained by Pinker and Laszlo (1992a). An effort to verify the 
above partitioning (Whitlock and LeCroy, unpublished results) showed the 65/35 ratio to be a good 
estimate for average atmospheric conditions. The PAR at the surface, F SPAR , was computed as 

Fspar = F sd | 0.42 + 2 ( zr □ 0.5 ) 2 } , (49) 

for both clear and cloudy conditions. It is believed that Eq. (49) was developed by fitting curves to the 
results derived by Pinker and Laszlo (1992b). 

4. Input data sources 

The most extensive prior application of this algorithm was for the 8-year period (July 1983 to 
June 1991) for which monthly average global SRB fields have been published (Gupta et al. 1999). All 
required cloud parameters, column precipitable water (PW), and column ozone for that work were taken 
from the ISCCP-C1 datasets. The latter two, namely, PW and column ozone were TOVS products 
incorporated into ISCCP-C1 datasets. Surface albedos were derived from literature formulas and from 
ERBE clear-sky TOA albedos, as discussed in Sec. 3.4. Aerosol properties used were climatological 
average values for four standard aerosol types, namely, maritime, continental, desert, and snow/ice 
(Deepak and Gerber 1983). The surface was classified as one of five types, namely, ocean, coast, land, 
desert, and snow/ice. A single aerosol type, or a combination of two, was associated with each of the 
surface types as described in Appendix C. Flux computations were made on a grid-box basis for the 
6596-box equal-area grid which has a resolution of about 280 km x 280 km. 

The state-of-the-art for some of the above datasets has advanced considerably over the last few 
years, and newer datasets are being used for the current applications. Cloud properties used in the current 
work are derived from pixel level ISCCP data, known as the DX data (Stackhouse et al. 2001) with the 
same algorithms as used for deriving the ISCCP-D products. Further, these cloud properties are being 
derived on an equal-area global grid, which consists of 44016 boxes and is called the nested grid. Areas 
of the boxes of this grid approximately equal the area of a 1 0 x 1 0 box at the equator. The increased 
spatial resolution provides insights into the structure of clouds not achievable with Cl or D1 datasets. 
Column PW for the current work is being derived from the data assimilation model products of the 
Goddard Earth Observing System, version- 1 (GEOS-1; Schubert et al. 1993), produced by the Data 
Assimilation Office (DAO) at the NASA Goddard Space Flight Center (GSFC). The GEOS-1 column 
PW is produced 6-hourly, and thus contains a representation of the diurnal variability. By contrast, the 
ISCCP-C1 column PW was a once/day product from TOVS, and contained no diurnal variability. 
Column ozone used for the current work comes from a long record available from the Total Ozone 
Mapping Spectrometer (TOMS), which flew aboard Nimbus-7 and Meteor-3, and is presently flying 
aboard EP-TOMS. The TOMS ozone product is deemed to be considerably superior to the ISCCP-C1 
column ozone, which is a TOVS product derived from an infrared channel on the HIRS-2 instruments. 
Column PW from GEOS-1 and TOMS ozone were both regridded to the nested grid to ensure 
compatibility with the new cloud products. Surface albedos and aerosol properties are still obtained from 
the same sources as in the earlier work. A value of 1365 Wm 2 for the solar constant, based on ERBE 
measurements was used hi the earlier work, and is also being used for the current work. 
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Table 1. Inputs required for LPSA and data sources used in the past, for the present work, and expected to be used 
in the future. 


LPSA input 

Past 

Data sources 
Present 

Future 

Solar constant 

ERBE 

ERBE 

TBD 

TOA reflected radiances 

ISCCP-C1 

ISCCP - DX 

ISCCP - DX 

Cloud amount 

ISCCP-C1 

ISCCP - DX 

ISCCP - DX 

Cloud optical depth 

ISCCP-C1 

ISCCP - DX 

ISCCP - DX 

Column PW 

TOVS 

GEOS - 1 

GEOS - 3 

Surface albedo 

LF, ERBE 

LF, ERBE 

TBD 

Column ozone 

TOVS 

TOMS 

TOMS 

Aerosol properties 

D&G 

D&G 

TBD 


Key: TBD - To be determined; LF — Literature formulas; D&G — Deepak and Gerber (1983). 

Datasets of still better quality are continuously coming online and will be used as they become 
available. For example, column PW for future work is likely to come from GEOS-3 or later versions of 
the data assimilation model now in use at the DAO. Newer models of aerosol spatial and temporal 
distribution and their optical properties are being explored for use in future work, as are the newer sources 
of surface albedo. TOMS is likely to continue as the future source of column ozone. Values of solar 
constant obtained from newer measurements will be examined to ascertain if their use in place of the 
ERBE-based value is warranted. A concise summary of the past, present, and future input data sources for 
LPSA is presented in Table 1. 

5. Results and discussion 

The current version of LPSA and the input datasets described in Sec. 4 have been used to derive a 
number of surface SW parameters for all months of 1986 and 1992. A complete list of these parameters 
is given below: 

1 . Clear-sky insolation, 

2. All-sky insolation, 

3. All-sky net SW flux, 

4. Direct SW flux, 

5. Diffuse SW flux, 

6. Photosynthetically active radiation, and 

7. All-sky surface albedo. 

Since the primary purpose of this report is to explain and document the scientific basis of the 
algorithm as much as possible, only small samples of these results will be presented and discussed here. 
Detailed presentations and discussions may be undertaken in the future when such results are produced in 
the context of various research projects. It suffices here to show that the results are physically consistent 
with the input meteorological fields used and with similar results from earlier work. With those 
objectives in mind, all-sky insolation results for 1992 are highlighted in this report. All-sky insolation is 
the most widely measured and used surface SW parameter. Also, the 1992 results are completely new for 
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this algorithm, being outside the July 1983 — June 1991 period, for which similar results were derived 
using ISCCP-C1 data for inputs (Gupta et al. 1999). The results for 1992 are compared with results for 
1986, both derived using the current version of the algorithm. This comparison may provide an estimate 
of interannual differences between 1986 and 1992, if any, as both sets are derived using identical input 
sources. Also, current 1986 results are compared with corresponding results derived earlier using the 
same algorithm but with ISCCP-C1 inputs. The latter comparison provides an estimate of the differences 
arising from (i) DX vs. Cl cloud inputs, (ii) GEOS-1 vs. TOVS meteorological inputs, and (iii) TOMS vs. 
TOVS ozone. 

Table 2. Comparison of hemispheric and global average all-sky surface insolation (Wm'“) for January, July, and 
the whole year for 1992 and 1986 from current work, and for 1986 from earlier work with ISCCP-C1 
inputs. 



N. H. 

S. H. 

Global 



1992 — Current Work 

Jan. 

128.5 

256.0 

192.3 

Jul. 

245.2 

116.6 

180.9 

Ann. 

191.2 

185.2 

188.2 



1986 - Current Work 

Jan. 

127.2 

246.7 

186.9 

Jul. 

243.6 

115.4 

179.5 

Ann. 

189.8 

184.4 

187.1 



1986- From ISCCP-C1 

Jan. 

126.8 

249.7 

188.2 

Jul. 

241.7 

118.0 

179.8 

Ann. 

186.6 

182.9 

184.8 


Figure 1 shows the seasonal variability of all-sky surface insolation averaged over the 
hemispheres and the globe for 1992 and the two datasets for 1986. All plots show a strong seasonal 
variability for the hemispheric averages, and a weak one for the global average. Table 2 presents numbers 
based on the same datasets for hemispheric and global averages for January, July, and the whole year. 
Seasonal variability (January to July difference) for the SH shows a slightly larger magnitude than for the 
NH. This difference arises from two reinforcing causes. First, the Sun-Earth distance is minimum during 
January (SH summer) and maximum during July, providing a stronger seasonal contrast over the SH. 
Second, there is a large seasonal variability of column water vapor in the NH, with a strong maximum in 
July (NH summer), which lowers the seasonal contrast for the NH. The numbers in Table 2 seem to 
indicate interannual differences between the 1992 and 1986 results from the current work, and those 
arising from the use of different input sources between the two datasets for 1986. It is emphasized that 
the above differences are presented only for illustrative purposes, and are not meant to establish 
interannual variability of surface insolation or characteristics of the input data sources. 
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Figure 2. Geographical distribution of monthly average surface insolation (W m 2 ) for 
January and July 1992 derived with inputs used in the current work. 
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Figure 3. Zonally-averaged surface insolation for January' and July for (a) 1992 from current work, 
(b) 1986 from current work, and (c) 1986 derived for ISCCP-C 1 inputs. 
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Figure 2 shows the geographical distribution of monthly average surface insolation for January 
and July 1992. The highest values occur over subtropical subsidence regions and polar areas of the 
summer hemisphere. Small cloud amounts over subsidence regions and longer sunshine duration over 
polar areas account for these features. Large cloud amounts in the stonn tracks account for the much 
lower values over midlatitudes of the summer hemisphere. For the winter hemisphere, polar areas do not 
receive sunlight, and the low values over midlatitudes result from a combination of low TOA insolation 
and dense cloud cover in the storm tracks. Relatively low values over the equatorial regions are caused 
by the heavy cloudiness along the intertropical convergence zone (ITCZ). 

Figure 3(a) shows zonally-averaged surface insolation for January and July derived from the 
same 1992 dataset. Figures 3(b) & 3(c) show corresponding results from the two datasets for 1986. The 
only significant differences between the three panels are seen in the January curves poleward of 60° in the 
SH. Qualitative comparisons of the input parameters used for computing these insolation fields showed 
that these differences arise primarily from corresponding differences in the cloud fields. Cloud amounts 
between 60° S and the South Pole in the 1986 dataset used in the current work were significantly higher 
than for the other two datasets. The results presented above demonstrate that the fields of SW parameters 
derived in the current work are physically understandable and consistent with the results of earlier work. 

A high degree of confidence in the quality of all-sky surface insolation shown above 
notwithstanding, a note of caution is in order regarding the quality of the direct, diffuse, and PAR fluxes 
derived with this algorithm. The fonnulas used to derive these products appear highly empirical and the 
results have never been validated. For these reasons, the direct, diffuse, and PAR fluxes were never 
distributed to the science community in the past, and will not be distributed in the near future. The 
inclusion of Eqs. (47) - (49) in this report was driven by the desire to present all elements of the 
algorithm as proposed by Staylor. The authors do not recommend the use of these formulas by the 
readers without independent validation. 

6. Planned enhancements 

A number of steps are being undertaken to update and enhance this algorithm to prepare it for 
newer applications. The first of these will be to change the time resolution from daily to instantaneous so 
that the results can be compared with high temporal resolution surface observations. The next will be to 
examine each of the attenuation terms in Eq. (23), including the models and measurements on which 
those terms are based. New parameterizations of the attenuation processes are currently under 
development by the authors and will be incorporated into this algorithm after they are validated. Newer 
models of aerosol spatial and temporal distribution and their optical properties are being explored as are 
newer sources of surface albedo. Those too will be incorporated into the algorithm when available. It is 
also planned to undertake a critical examination of Eqs. (47) — (49) which are used for computing direct, 
diffuse, and PAR fluxes respectively. Results of these equations will be compared with other model 
results and observations, and the equations will be modified as necessary. 
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Appendix A 

Astronomical Relationships 


The insolation reaching the Earth is governed by the inverse square law through the eccentricity 
correction factor, (djdf, in Eq. (2). The mean Sun-Earth distance, d m , is 1.496 x 10 11 m and is called an 
astronomical unit (AU). The instantaneous Sun-Earth distance varies from 1.471 x 10 11 m (0.983 ALT) in 
early January to 1.521 x 10 11 m (1.017 AU) in early July. Iqbal (1983) presents several simple 
expressions for the daily average value of (djdf and recommends one derived by Spencer (1971): 

( djd f = 1.000110+ 0.034221 cosT + 0.001280 sin T (A1) 

+ 0.0007 19 cos 2T+ 0.000077 sin IT, 


where T (in radians) is given by 


r = 20(4, Gl)/ 365, (A2) 

for a year of 365 days. Equations (Al) and (A2) were used in the present work with 366 substituted in 
Eq. (A2) for the leap years. According to Iqbal (1983), the maximum error in ( djdf computed with Eq. 
(Al) is 0.0001. 

Another astronomical variable which affects insolation reaching the Earth is the solar declination, 
0, through cos Z in Eq. (3). Solar declination varies from +23.5° at the summer solstice (about 21 June) 
to -23.5° at the winter solstice (about 22 December) and goes through zero at the vernal and autumnal 
equinoxes (about 21 March and 22 September, respectively). Note that the above seasonal descriptions 
apply to the NH; the opposite apply to the SH. Iqbal (1983) presents several expressions for the daily 
average value of 0and recommends one, also derived by Spencer (1971): 

0 = ( 0.006918 □ 0.399912 cosT + 0.070257 sinT 

□ 0.006758 cos 2 T + 0.000907 sin IT (A3) 

□ 0.002697 cos 3T + 0.00 1 48 sin 3T ) ( 1 80 / 0 ) , 

where T is already defined in Eq. (A2). Equation (A3) was used in the present work. According to Iqbal 
(1983), the maximum error in 0 computed with Eq. (A3) is 0.0006 radian. 
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Attenuation Factor 


Appendix B 

Attenuation Under Clear Skies 

The absorption/scattering processes of the various constituents (excluding aerosols) contributing 
to the attenuation of solar radiation in the atmosphere are represented by Eqs. (24) — (28) for an overhead 
Sun. As stated in Sec. 3.2, equations in these exact forms were not found in the cited references, namely, 
Lacis and Hansen (1974), and Yamamoto (1962). However, functionally equivalent equations 
representing attenuation due to water vapor, ozone, and Rayleigh scattering are given in the Lacis and 
Hansen reference. Comparisons of Eqs. (24) and (25) with corresponding formulas from Lacis and 
Hansen are presented below in Figs. B1 and B2 respectively. Figure B1 shows good agreement for 
column water vapor values above 0. 1 pr-cm which makes the use of Eq. (24) appropriate over most of the 
globe. Figure B2 shows good agreement thoughout the range. 

Absorption due to C0 2 and 0 2 (Eqs. 26 and 27) is related to surface pressure, P s . At sea level (P s 
= 1), attenuation by C0 2 amounts to about 0.63% of the TOA insolation which is close to that obtained 
from the curves for C0 2 absorption given by Manabe and Strickler (1964). Absorption by 0 2 amounts to 
about 0.75%. This value agrees well with those obtained from: i) the Hoyt (1978) formula, ii) the 
parameterization by Chou and Suarez (1999), and iii) line-by-line calculations made by the authors 
(unpublished results). Attenuation due to molecular scattering (Eq. 28) amounts to about 3.5% of the 
TOA insolation. This number is in close agreement with that obtained from Eq. (41) of Lacis and Hansen 
which represents atmospheric albedo due to Rayleigh scattering. Also, it is important to note that Eq. (28) 
represents only the upward part of the scattered radiation. The downward part remains in the downward 
radiation stream. 
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Figure B 1 . Attenuation by water vapor. 


Figure B2. Attenuation by ozone 
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Appendix C 

Aerosol Distribution and Radiative Properties 

The magnitude of the aerosol attenuation in this model was computed on the basis of five surface 
scene types listed in Table Cl. Standard aerosol types (e.g., maritime, continental, desert) and values of 
do, and g for them were adopted from Deepak and Gerber (1983) with minor adjustments. A standard 
aerosol type (or a combination of them) was associated with each surface scene type. In addition, 
climatological mean values of den further parameterized by Staylor in terms of u for each surface scene 
type, were used. All of this information is presented in Table Cl below. 

Table Cl. Aerosol radiative parameters fl,,,., /7 0 , and g for the five surface scene types. A Tch . in the formula for £J fr is 

the ERBE clear-sky TOA albedo. 


Scene type 

Aerosol type 

ELer 

Do 

g 

Ocean 

Maritime 

0.15 m 

0.98 

0.60 

Land 

Continental 

0.35 m 

0.90 

0.66 

Desert 

Desert 

(03+0.5A Tclr )u 

0.92 

0.60 

Coast 

50/50 Maritime & 

0.25 u 

0.94 

0.64 


Continental 




Snow/Ice 

Snow/Ice 

0.03 

0.97 

0.67 
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Appendix D 

Sampling Procedures and Alternative T c Computations 

The regression coefficients Cj and C 2 (Eq. 35) for computing overcast reflectance, R orc , for each 
grid box and each day were determined by linear regression between R oyc fx^i 0 and \u li 0 I(u+ijl 0 )\’ . 
Regression was performed separately for each ISCCP satellite for every month. Only overcast grid boxes 
for which solar and view zenith angles were less than 70° were sampled for regression. This procedure 
still left datasets of more than 100,000 grid boxes, a size that was considered too cumbersome for 
regression analysis. These were reduced to a more manageable size (less than 10,000) by sorting R mCl uu„ 
and on cloud optical depth and retaining only those in the highest 10 — 20% of optical depth 

range. 

Corresponding coefficients for computing clear-sky reflectance, R dn over ice-free oceans (C 3 and 
C 4 in Eq. 36) were determined by linear regression between R dr and (uu„) °' 75 . Clear grid boxes, subject to 
the same restrictions of solar and view zenith angles as above, were sampled for this regression. Many of 
these datasets also had more than 40,000 values, still too large for regression analysis. These datasets 
were randomly sampled to reduce the number down to about 4,000. 

Two conditions encountered under which the use of Eq. (34) was deemed inappropriate for T c 
computation were the following: 

1) When (R orc - R dr ) was too small. A lower limit of 0. 15 was imposed on (R mc - R dr ). 

2) When (R mc - R mea J was negative. Under these conditions, when both cloud amount and cloud optical 
depth were available, T c was computed as 

T c = 0.05 + 0.95 (1[] 0.2 A c [] c 031 ), ( D1 ) 

where A c is the fractional cloud amount and jj : is the cloud optical depth. When only A c was available, T c 
was computed as 


T c = 0.2 + 0.8 (1 G A c ) 07 . (D2) 

Both equations provide T c = 1 for A c = 0 (clear sky). For A c = 1 (overcast), Eq. (Dl) provides the 
minimum value of T c (= 0.05) for a jj : value of about 80. This was assumed to be the highest average 
value [J: for a dense overcast. Also for A c = 1, Eq. (D2) provides T c = 0.2, which in terms of Eq. (Dl) 
corresponds to a value of /J of about 50. 
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Appendix E 

Alternative Surface Albedo Computations 

As stated in Sec. 3.4, when snow/ice cover was present in a grid box, alternative methods were 
used for computing A Sclr and A Sm . c as follows: 

Oceans: When an ERBE-derived value of A Sclr was available, it was used as such and a 
corresponding value of A Sm . c was computed as 

A™ = 4*(“/ °-6) (1Di) , (El) 

where s represents the fractional snow/ice cover for the grid box. When A Sclr from ERBE was not 
available, it was computed as 


A Sclr = A,(lDs) + 0.5s, (E2) 

where A Sl represents a value computed from Eq. (39) but capped at 0.25. The corresponding value of A Soyc 
was computed as 


4* = 0.065 (1 0^) + 0.5 5. (E3) 

Eq. (E2) reduces to Eq. (39) when 5 = 0 (ice-free ocean) but caps the value of A Sch . at 0.25, and Eq. (E3) 
reduces to Eq. (40). Both (E2) and (E3) cap surface albedo values at 0.50 when 5=1. 

Other surfaces: When snow/ice cover is present in a grid box and an ERBE-derived value of 
surface albedo is available, that value is used for both A Sdr and A Sovc . When an ERBE-derived value is not 
available, A Sdr is computed as 

A Scll . = 0.2 (1 [] 5 ) + 0.7 5 , (E4) 


and the same value is used for A Sovc . 
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